NASA Technical Memorandum 101298 


Discretization Formulas for 
Unstructured Grids 


(USA-TK-io 12S8 


tlSCfii 1 1ZA1 1C 
CHICS ( hASt) 


h 


iCBaULAS 
12 pcs Cl 20B 


was - 10242 


. Unclas 
g3 /34 0166358 

Kenneth J. Baumeister 
L^vm Research Center 
Cleveland, Ohio 


Prepared for the 

Second International Conference on Numerical Grid Generation 
in Computational Fluid Dynamics 

cosponsored by the National Aeronautics and Space Administration, 

the U.S. Air Force Office of Scientific Research, and the University of Miami 

Miami Beach, Florida, December 5-8, 1988 





DISCRETIZATION FORMULAS FOR UNSTRUCTURED GRIDS 


Kenneth J. Baumeister 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


ABSTRACT 

The Galerkin weighted residual technique using linear tri- 
angular weight functions is employed to develop finite differ- 
ence formula in cartesian coordinates for the Laplacian 
operator, first derivative operators and the function for 
unstructured triangular grids. The weighted residual coeffici- 
ents associated with the weak formulation of the Laplacian 
operator are shown to agree with the Taylor series approach on 
a global average. In addition, a simple algorithm is presented 
to determine the Voronoi (finite difference) area for an 
unstructured grid. 

INTRODUCTION 

In developing unstructured finite difference equations, 
Jameson [1, Eq. (4.10)] has applied the standard weighted 
residual Galerkin method to obtain a time dependent discretiza- 
tion of the Euler equations for an unstructured triangular 
mesh. Erlebacher [2] has utilized variational methods on a 
pointwise basis to establish a difference operator for the 
Laplacian operator for a central difference cell similar to 
that shown in Fig. 1. The present paper will utilize the 
global Galerkin weighted residual technique to also generate 
simple closed form finite difference approximations for the 
Laplacian operator, first derivative operators, and the func- 
tion itself for the general central cell in Fig. 1, as well as 
the boundary element cell also shown in Fig. 1. The boundary 
element cells are required in acoustics and electromagnetic 
theory in establishing the far field radiation conditions. 

GALERKIN WEIGHTED RESIDUAL FORMULATION 

Consider the following partial derivatives: 



V • I^V* ; u 2 ff ; & 3 If ; V (1) 

4» is a scaler or potential quantity and the |3's are variable 
property coefficients. To obtain the finite difference expres- 
sions for the terms in Eq. (1), the continuous domain D is 
first divided into discrete triangular areas A e staked out by 
nodal (grid) points P as shown in Fig. 1. The number of con- 
nected areas (called a cell) needed to define the difference 
equation for node o is labeled M. 



The difference equations developed in Appendices A and B 
take on the form 
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where the Laplacian operator is converted to the later form by 
means of the weak formulation [3, p. 4433 of the weighted 
residual approach. The values of the coefficients k and a 
depend only on the location of the grid points and the proper- 
ties 3- The superscript * represents either the derivative 
in x ( * = x) or y (* = y) and the variable itself (* = 4>). 

The coefficients k* for $ and its first derivatives are 
included in Appendix C. 
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The major effort of the present work is to establish the 
relationship between the Laplacian operator and the a coeffi- 
cients from the method of weighted residuals. An expression 
for the Laplacian of the form 

y2<j> = f (ai , a2 , • • • ) (4) 

is desired. Using the result from Appendix D for a linear 
shape function, the Laplacian can be expressed as 

2 Vo + Vl + a 2*2 + ••• + Vo 

V^<j> = -2_o \ — ! — E_E + r (5 ) 
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where Ay is the total area of the cell. The a values in 
Eq. (5) are given as 


% - -[* <nm, '«A n J[(v - 0 

( y o - v) + ( x n - v)( x n- - x o). 
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where the subscripts are defined in Appendix C. The term r 
is the contribution from the natural boundary condition that 
comes automatically from the weak formulation of the method of 
weighted residuals and contributes only on the boundary of the 
doma i n . 

DISCUSSION OF RESULTS 
Uniform Grids 

For the conventional six-point hexagon difference grid as 
shown in Fig. 2, the difference equations for the Laplacian, 
first derivatives, and the function as calculated from Eqs. (2) 
and (5) are in agreement with the standard difference operators 
that appear in the literature [4, p. 1114]. Similarly, the 



very popular six-node system shown in Fig. 3 also shows agree- 
ment. For brevity only the Laplacian operator will be shown in 
Fig. 3 and on the remaining figures. 


As is well known in the literature [3, p. 105], the method 
of weighted residual and the conventional Taylor series expan- 
sion often yield different results. For the conventional four- 
point square difference grid as shown in part A of Fig. 4, the 
difference equation for the Laplacian is different by a factor 
of 3/2 from the conventional Taylor series representation of 
the Laplacian operator [4, p. 1114]. Similarly, the eight-node 
system shown in part B of Fig. 4 differs by a factor of 3/4. 
However, since the method of weighted residual (finite element 
theory) is guaranteed to converge, these are acceptable formu- 
lae provided they are applied to a consistent global mesh and 
not a single arbitrary cell. 

Relationship Between FE and FD Methods 


For uniform grids, of the type displayed in Figs. 2 to 4, 
the finite difference and finite element expressions for the 
Laplacian operator will be in agreement with the conventional 
Taylor series difference equation provided that the bias con- 
stant B, defined by the following area rule, is identical to 
zero. 
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FIGURE 2. - UNIFORM 6 MODE HEXAGON GRID SYSTEM. 



FIGURE 3. - UNIFORM 6 NODE SLANTED GRID SYSTEM. 
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FIGURE - GLOBAL CONVERGENCE EXAMPLE WITH RECTANGULAR ELEMENTS. 


B = 3A V /A T - 1 (9) 

A y is the Voronoi neighborhood [4] or the area normally 
employed in the finite difference analysis. If B ^ 0, the 
coefficient in front of the finite element expression will dif- 
fer from the expression for the Taylor series expression by a 
factor of 1 + B. If B = 0, a single pattern of uniform ele- 
ments can normally be placed in a regular domain. 

Global Relation To Taylor Series 

Although the difference equations for the Laplacian opera- 
tor in Fig. 4 differ in each cell from the standard Taylor 
series difference equation, on a global average the Taylor 
series and finite element equation will average out to the same 
value. In this case the area average of the eight-node ele- 
ment combined with the four-node element will be such that the 
Taylor series will be valid. As shown in Fig. 4, the average 
global value of B is equal to zero. (B is a measure of the 
difference between the finite element and Taylor series as just 
defined.) A similar situation occurs for triangular grids, 
although the results are not shown. 

Nonuniform Grids 

Figure 5 displays the difference equations for a nonuni- 
form five-node cell. The Laplacian operator shown in Fig. 5 
is only valid when taken in conjunction with the other grid 
systems that surround it in a given domain. It is not to be 
thought of as a Taylor series approximation for the cell shown 
In Fig. 5. 
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FIGURE 5. - NONUNIFORM 6 NODE GRID SYSTEM. 

Voronoi (Finite Difference) Area 

The determination of the Voronoi area plays a major role 
in the mesh generation by triangulation of arbitrary points 
[2]. The coefficients used to generate the Laplacian operator 
can also be used to easily determine the Voronoi area without 
geometric construction; 

. «o( x o * 4) *]( x i * 4) * 2(4 ♦ 4) 

A V 4 + 4 + 4 + --‘ 



where a's are given by Eqs. (6) and (7). The rationale for 
development of this simple algorithm to predict the Voronoi 
area comes from the fact that replacing Ay/3 by Ay in 
Eq. (5) fortuitously satisfies the Taylor series approximation 
for uniform grids. (However, it does not in general for non- 
uniform grids, and Ay/3 = Ay when B = 0.) Therefore, if $ 
is assumed to vary as [x 2 + y 2 ]/4 then the Laplacian for uni- 
form grids will have a value of unity and Eq. (10) results. 

To validate this hypothesis, Eq. (10) was checked against the 
nonuniform cell shown in Fig. 6, as well as a large variety of 
other nonuniform cells, and found to always be in exact agree- 
ment with geometrical calculations. For certain complex cells, 
the Voronoi boundaries can cross as shown in Fig. 6. In these 
cases, the algorithm will consider reversed areas as negative, 
in that it senses the direction of the path of the Voronoi 
region. Again the algorithm agrees with the geometric 
calculation. 
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AXIAL COORDINATE, X 

FIGURE 6. - CALCULATION OF VORONOI AREA FOR 
G NODE GRID WITH REVERSAL. 


Local Convergence 

If B = 0, Eq. (5) combined with Eq. (9) for central cells 
yields 

= (ao4>o + a ]$i + a 2^2 + ••• + a p4>p)/A\/ (ID 

where Ay is the Voronoi area of the cell and the a values 
are given in Eqs. (6) and (7). Equation (11) was checked 
against a number of nonuniform grids (six or less nodes) where 
B * 0, and found to satisfy the test functions 1, x, y, x 2 , 
and y 2 . However, the test function xy was not satisfied 
when B * 0. For the special case, as shown in Fig. 5, where 
B ~ 0, Eq. (11) nearly satisfied the xy function test. 
Therefore, Eq. (11) may be valid on a pointwise bases for the 
special condition of B = 0. More testing will be required 
for val idation. 

CONCLUSIONS 

The Galerkin weighted residual technique using linear 
triangular weight functions is employed to develop finite dif- 
ference formulae in cartesian coordinates for the Laplacian 
operator, first derivative operators and the function for 
unstructured triangular grids. 
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APPENDIX A - DERIVATION GALERKIN WEIGHTED RESIDUAL EQUATIONS 

In accordance with the method of weighted residuals, 
weight functions W 0 in the form of a linear pyramid are 
Introduced (one for each grid point), such that the Integra- 
tion of the product of each term in Eq. (1) times the weight 
W 0 over the cell area shown in Fig. 1 yields 


:-VW 1 • (3 1 V<t.)dA = ^ 


,(e) 


e=l 


-VN< e 2[v (e >] 4 (e) } 


dA 


(Al) 


W 1 P 2 ff dA = X P 


(e) 

2 


e-i 




dA (A2) 


where the subscript x stands for the x derivatives of the 
known shape functions and N stands for .the linear interpola- 
tion functions whose x and y dependence is presented in 
most finite element texts [3, p. Ill, Eq. (3.1.4)]. The equa- 
tions for the derivative of <j> with respect to y and the 
function itself are similar to Eq. (A2). The Laplacian opera- 
tor in Eq. (Al) has utilized the weak formulation of the 
weighted residual formulation C 3 , p. 443]. 

APPENDIX B - AVERAGED OPERATORS 


Assuming averaged finite difference values over the 
weighting function area allows the function and first deriva- 
tives to be written explicitly. For example. 


a 2 tf"“ 15 

' cl 


(e) 

2 




e). 


dA 


(Bl) 


In this case, the integral of the weight W 0 over the 
cell area is equal to Aj/3, which was determined with the aid 
of the standard area integration formula C3, p. 112]. Equa- 
tion (Bl) as well as the equation for the Laplacian, the func- 
tion and the derivative of the function with respect to y 
can now be easily evaluated using conventional finite element 
theory to obtain the difference equations given in Appendix C 
and Eqs. (6) and (7) in the body of the report. 

APPENDIX C - COEFFICIENTS FOR FUNCTION AND FIRST DERIVATIVES 

To develop the difference equations, the Kronecker delta 
is defined as 

^a,p = 0, a * p; or &a,p = 1 » a = p (Cl) 
and the following circular indices are also defined 

n- ■ n - 1 + PS , (C2) 

n , 1 

n+ = n + 1 - P$ n p (C3) 

e+ = e + 1 - D (C4) 

e.P 

nm = n - 1 + M6„ , (C5) 

n, l 


wi th 


A e = [1/2] x 0 (yj - y k ) + [1/2] xj(y k - y 0 ) + [1/2] x k (y 0 - yj> 

(C6> 

[3, p. 110, Eq. (3.1.2)]. Using this notation, the coeffici- 
ents for the following operators can be written as follows: 

x derivative operators 



+ 



T+l-n ,P 


,) 
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(C8) 
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derivative operators 


k n ' [^3° ’' 2 A t](v x o)( 1 ‘ 5 d,P-m) 

♦ (S‘ n> /2A T )(K 0 - ,J[t - ‘p*1-„.p_ M ] (C9> 

k 0 ' g ( B 3 3>/2A T)[V- %] < C,0 > 

4> function 

< - K nm> ( A n m /4A T)](' - £ n,P-n) 

* K"’ (V 4 A t)](' - *P+1 -n . P- m) <Cm 

M 

k .♦ = (1/2 A t ) ^ (3^ e) A e (Cl 2) 

e=l 

APPENDIX D - SOLUTION OF THE WEIGHTED RESIDUAL EQUATIONS 

Consider the partial differential equation of the form 

V 2 4> + !£- = 0 (Dl) 

Employing the weak solution [3, p. 443] for the four-node cell 
in part A of Fig. 4 yields 

(4>1 + 4>2 + 4>3 + $4 - 4<j> 0 ) + h [ ( 4> i - 4>3 ^ /3 ] = 0 (D2) 

To put Eq. (D2) in a more familiar form which coincides with 
the standard Taylor series finite difference approach, Eq. (D2) 
is multiplied by 3/Aj, which corresponds to 3/2h 2 for the 
geometry shown in Fig. 4. Thus, Eq. (D2) becomes 

3/2 [<4>i + $2 + $3 + 4*4 - 44> 0 )/h 2 3 + (<J>i - 4>3 ^ /2h = 0 (D3) 

The 3/Ay term was chosen as the multiplying constant because 
it will always convert the function and its first derivatives 
so that they match the Taylor series expression. Consequently, 
the Ay/3 term appears in Eq. (5) in the body of this report. 
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